function Plot_Solution_Model_20220926(parms,gesol,prob_Nxz)

    if nargin==3
        bUseDefaultRange = false;
        % select relevant range
        P1 = 0.001; P2 = 0.999;
        idx_N1 = Inf; idx_N2 = -Inf;
        idx_x1 = Inf; idx_x2 = -Inf;

        for iz = 1:parms.Nz
            cdf_N = cumsum(sum(prob_Nxz(:,:,iz),2)); cdf_N = cdf_N/cdf_N(end);
            idx_N1 = min(idx_N1, find(cdf_N<=P1,1,'last'));
            idx_N2 = max(idx_N2, find(cdf_N>=P2,1,'first'));
            cdf_x = cumsum(sum(prob_Nxz(:,:,iz),1)); cdf_x = cdf_x/cdf_x(end);
            idx_x1 = min(idx_x1, find(cdf_x<=P1,1,'last'));
            idx_x2 = max(idx_x2, find(cdf_x>=P2,1,'first'));
        end
    else
        bUseDefaultRange = true;
    end

    m = 2; n = 2; idxplot = 0;

    phi_of_x = @(x)1./(1+exp(-x));
    [NN,PP] = meshgrid(parms.grid_N,phi_of_x(parms.grid_x));
    
    % plot D/P ratio
%     DP = gesol.C./gesol.stock.price_cons_unlev;
%     gesol.DPratio = DP;
    
    % plot_list = {'S','CV','theta','C','V','f','g','wages','DPratio'};
    
    plot_list = {'S','CV','theta','C','V','f','g','wages'};
    
    for i = 1:numel(plot_list)
        temp = gesol.(plot_list{i});
        zlim1 = floor(min(temp(idx_N1:idx_N2,idx_x1:idx_x2,:),[],'all')*10)/10;
        zlim2 = ceil(max(temp(idx_N1:idx_N2,idx_x1:idx_x2,:),[],'all')*10)/10;
        idxplot = idxplot + 1;
        if idxplot>m*n; idxplot=1; end
        if idxplot==1; figure; end
        subplot(m,n,idxplot);
        mesh(NN,PP,gesol.(plot_list{i})(:,:,1)')
        xlabel('N')
        ylabel('\phi')
        title([plot_list{i},', z = 1'])
        set(gca,'zlim',[zlim1 zlim2])
        if ~bUseDefaultRange
            set(gca,'xlim',parms.grid_N([idx_N1 idx_N2]));
            set(gca,'ylim',phi_of_x(parms.grid_x([idx_x1 idx_x2])));
        end
        
        idxplot = idxplot + 1; subplot(m,n,idxplot);
        mesh(NN,PP,gesol.(plot_list{i})(:,:,2)')
        xlabel('N')
        ylabel('\phi')
        title([plot_list{i},', z = 2'])
        if ~bUseDefaultRange
            set(gca,'xlim',parms.grid_N([idx_N1 idx_N2]));
            set(gca,'ylim',phi_of_x(parms.grid_x([idx_x1 idx_x2])));
        end
        set(gca,'zlim',[zlim1 zlim2])
    end
    
    % plot zeta
    figure
    subplot(1,2,1)
    zeta_fun = @(N,phi)(N + (1-N).*(phi.^(1-1/parms.chi))).*((N + phi.*(1-N)).^(parms.gamma-1));
    zeta_NP = zeta_fun(NN,PP);
    mesh(NN,PP,zeta_NP)
    xlabel('N')
    ylabel('\phi')
    title('\zeta(N,\phi)')
    if ~bUseDefaultRange
        set(gca,'xlim',parms.grid_N([idx_N1 idx_N2]));
        set(gca,'ylim',phi_of_x(parms.grid_x([idx_x1 idx_x2])));
    end
    set(gca,'zlim',[floor(min(zeta_NP,[],'all')) ceil(max(zeta_NP,[],'all'))])
    subplot(1,2,2)
    mesh(NN,PP,log(zeta_NP))
    xlabel('N')
    ylabel('\phi')
    title('log \zeta(N,\phi)')
    if ~bUseDefaultRange
        set(gca,'xlim',parms.grid_N([idx_N1 idx_N2]));
        set(gca,'ylim',phi_of_x(parms.grid_x([idx_x1 idx_x2])));
    end
    set(gca,'zlim',[floor(min(log(zeta_NP),[],'all')) ceil(max(log(zeta_NP),[],'all'))])
    
    % plot SDF
    if parms.Nz==2
        
        figure; m = 2; n = 2; idxplot = 0;
        
        if bUseDefaultRange
            zmin = 0.1*floor(min(gesol.spd1period(:))*10);
            zmax = 0.1*ceil(max(gesol.spd1period(:))*10);
        else
            zmin = 0.1*floor(min(gesol.spd1period(idx_N1:idx_N2,idx_x1:idx_x2,:),[],'all')*10);
            zmax = 0.1*ceil(max(gesol.spd1period(idx_N1:idx_N2,idx_x1:idx_x2,:),[],'all')*10);
        end
        
        for iz = 1:2
            for iznext = 1:2
                idxplot = idxplot + 1;
                subplot(m,n,idxplot)
                mesh(NN,PP,gesol.spd1period(:,:,iz,iznext)')
                xlabel('N')
                ylabel('\phi')
                title(['SDF(N,\phi,z,z''), z=',num2str(iz),', z''=',num2str(iznext)])
                set(gca,'zlim',[zmin zmax])
                if ~bUseDefaultRange
                    set(gca,'xlim',parms.grid_N([idx_N1 idx_N2]));
                    set(gca,'ylim',phi_of_x(parms.grid_x([idx_x1 idx_x2])));
                end
            end
        end
        
    end
    

end